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£h ! Abstract 

t^. ■ We describe a method of solving the nuclear Skyrme-Hartree-Fock problem by using 

a deformed Cartesian harmonic oscillator basis. The complete list of expressions required 
to calculate local densities, total energy, and self-consistent fields is presented, and an 
implementation of the self-consistent symmetries is discussed. Formulas to calculate ma- 
trix elements in the Cartesian harmonic oscillator basis are derived for the nuclear and 
O ' Coulomb interactions. 
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1 Introduction 



Self-consistent methods have been used in the low-energy nuclear structure studies over 
many years, and represent a mature field with numerous successful applications, see Refs. 
|T], |[] for a review. A number of computer codes solving the nuclear Hartree- 

Fock (HF) problem have already been developed. Two types of effective nucleon-nucleon 
^ ■ interactions have been mainly employed. Starting with the work of Vautherin and Brink 
H many authors have applied the nuclear HF theory with the Skyrme effective interaction, 
while the work of Gogny [|7|, § initiated numerous studies with the force which carries his 
name. 

The methods employed to solve the HF equations depend mainly on the effective force 
used and on the assumed symmetries of the many-body wave functions. For the solutions 
which allow at least triaxial deformations, two different methods have been applied for 
the two abovementioned effective interactions. The first one, used in conjuncture with 
the Skyrme interaction, is formulated in the spatial coordinates and makes use of the 
finite-difference M, or Fourier [I0~|, or spline-collocation [O, [I2J methods to approximate 



differential operators. The solution is then obtained by using the imaginary time evolution 



operator |L3 



The second one, used for the finite-range Gogny interaction, employs a truncated 



harmonic oscillator (HO) basis [14, 15] and solves the problem either by an iterative 
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diagonalization of the mean- field Hamiltonian, or by the gradient- [16] or the conjugate 
gradient- jl7| methods. A similar basis-expansion technique has also been recently used 
for studies in the frame of the relativistic mean field theory, for a review cf. Ref. flT9 



The present study aims at describing the method which incorporates the advantages of 
both existing approaches, and combines the robustness of the Cartesian HO basis with 
the simplicity of the Skyrme interaction. 

The methods using spatial coordinates have several advantages. First of all, various 
nuclear shapes can be easily treated on the same footing; the same cubic lattice of points 
in three spatial dimensions is suitable to accommodate wave functions with, in principle, 
arbitrary deformations restricted only by a specific symmetrization of the lattice. This 
allows easy studies of systems for which the deformation is not a priori known, or is ill 
defined because of deformation instabilities or a shape coexistence. Secondly, using spatial 
coordinates allows studies which address the asymptotic form of nucleonic wave functions 
at large distances. This is particularly important for a precise description of weakly bound 
nuclei, where the use of spatial coordinates is a necessity f20| . Third, for the Skyrme zero- 
range interaction, the mean fields are local (apart from a velocity dependence) and can 
be easily programmed in the spatial coordinates. Last but not least, the treatment of 
wave functions on large lattices (12x12x12 is a typical example) is easily amenable to 
vectorization or parallelization of the algorithm. 

Methods using the HO basis have other advantages. Firstly, the basis provides a 
natural cut-off for many operators which otherwise are unbound and require particularly 
delicate treatment in the spatial coordinates. This concerns in particular the multipole 
moment and the angular momentum operators which are often used as constraining op- 
erators. For the corresponding constraints the solutions can become unstable when the 
non-zero probability amplitudes (wave functions) move towards large distances as it is the 
case for e.g. weakly bound nucleons. Secondly, much smaller spaces are usually required 
to describe the nuclear wave functions within a given precision. Typically a basis of about 
300 HO wave functions is sufficient for most applications. Third, the iterative diagonal- 
ization of the mean field Hamiltonian can be used to find the self-consistent solutions, 
which provides a rapidly converging algorithm, and, last but not least, scalar or super 
scalar computers can also be used, because the typical sizes of the information handled is 
smaller and the performance is less dependent on the use of a vector processor. 

The above listed advantages of a given method correspond often (although not always) 
to the disadvantages of the other one, and both methods described above are in this respect 
complementary. Our motivation to construct the Skyrme-Hartree-Fock code using the 
Cartesian HO basis was based on the necessity to obtain a tool which would allow rapid 
solutions for the nuclear superdeformed or hyperdeformed rotating states for which the 
deformation is relatively well known a priori. Indeed, the method employed by us gives 
a particularly fast, stable, robust, and simple algorithm to solve such physical problems. 

The paper is organized in the following way. In Sec. ^] we present the HF method 
for the Skyrme interaction, and in particular we discuss the local densities, total energy, 
constraints, mean fields, and HF equations. In Sec. ^] we describe the use of various 
symmetries and in Sec. [| the use of the Cartesian HO basis. Sec. |^ is devoted to the new 
method we use to calculate the direct Coulomb potential. The method described in the 
present study is implemented in the computer code hfodd, published in the subsequent 
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paper [EI] which is referred to as II. 



2 Hartree-Fock method 

In this Section we present the complete set of formulas which are used when solving the 
Skyrme-Hartree-Fock problem. 

2.1 Local densities 

In the HF approximation, the total energy of a nuclear system is, in general, a functional 
of the one-body non-local density defined as 

p a (rcr,r'cr') = (^\al l(T , a a T(Ta \^) . (1) 

Here, |^) is a many-body wave function, and a\ aa and a raa are the operators creating 
and annihilating a neutron (a=n) or a proton (a=p) in the space-point r and having the 
projection of spin cr=±|. 

The non-local density (fl|) can be written |£2J as a sum of the scalar, p a (r,r'), and 
vector, s a (r,r') terms, defined by 

Pa(r, r') = Pa(ra, r'a) (2) 

and 

s a (r, r') = Pa(ra, r'a') < a'\cr\a > (3) 

as 

p a (r<7, r'a') = -p a (r, r')5 a(T > + -s a (r, r') ■ a aa >. (4) 

Assuming that the total energy depends only on local (r=r / ) densities, and on their 
derivatives at r=r' up to the second order, one has to consider the following real nucleonic 
densities fl2l 



particle-density p a {r) and spin density s a (r), 

Pa(r) = p a (r,r), (5a) 

s a (r) = s a (r,r), (5b) 

kinetic density r a (r) and vector kinetic density T a (r), 

r Q (r) = [V • Vpjr, r')] r=r , , (6a) 

T a (r) = [V- V'a a (r,r')] r=r , , (6b) 

momentum density j Q ,(r) and spin-current tensor J a (r), 

i„W = ^[(V- V')p Q (r,r')] r=r , (7a) 

^a(r-) = ^[(V M -V;) S , )a (r,r')l ^ . (7b) 



r=r' 
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In what follows we often omit the space argument r of local densities. 

For each of the local densities we define the isoscalar and isovector density. For exam- 
ple, the isoscalar (p ) and isovector (pi) particle densities are defined as 

PO = Pn + Pp, Pi = Pn - P V - (8) 

Following a traditional notation, for isoscalar densities we sometimes omit the isospin 
index, for example, p=po- 



2.2 Hartree-Fock energy 

In the Skyrme-HF approximation, the total energy £ of a nucleus is given as a sum of the 
kinetic, Skyrme, Coulomb, and pairing terms 

g £k' n _|_ ^Skyrme _|_ ^Coul _|_ £pair 

The kinetic energy of both protons and neutrons is given by the integral of the isoscalar 
kinetic density T (r), cf. Eqs. ( |oa| ) and (|j), 

^l^-i)/^' < 10 > 

where the standard factor (1— 4) provides a simple approximation to the center-of-mass 
correction [23[ in terms of the number of nucleons A. 



The Skyrme energy is the sum of the isoscalar (t=0) and isovector (t=l) terms, and 
is given as the integral of two energy densities. The first of these densities, 7if en (r), 
depends on the time-even densities p t , r t , and J t = J^u while the second one, 7i° dd (r), 



depends on the time-odd densities, s t} T t , j t (see Ref. |23|]), i.e. 



^Skyrme = J- J £r (ftf^r) + ft^V)) » ( U ) 



t=0,l 



for 



Ur n (r) = C?p 2 t +C^p t Apt + C; Pt T t + C t J J t +C7 J p t V-J t , (12a) 
Hf d (r) ee C°s 2 t +C t As s t -As t + C?s t .T t + Cij 2 t +C7 3 s t -(Vxj t ). (12b) 

2 

In Eq. (p.2a|), the square of the tensor density is defined as J t =J2nu Juvti an d its vector 



part J 4 is defined as J\,t=T,^ e x^uJ^u,t- 

Expressions introduced via Eqs. (12), relating the ten time-even coupling constants 
Cf, C t Ap , C t r , C/, and C t VJ 5 an d the ten time-odd coupling constants Cf, C t As , C t T , Cf, 



and C t \ to the standard parameters of the Skyrme interaction are given in Ref. [2~4"| 
Internally, the code hfodd uses the coupling constants in the traditional representation 
in which every term is a sum of the total density squared, and of the sum of squares of 
neutron and proton densities. Every such a term can also be written down in the isoscalar- 
isovector representation used in Eqs. (|TTD and (12). For the simplest term depending on 
particle densities, the correspondence is 

CL(Pn + P P f + C? um (pl + pi) = C p ( Pn + Pp f + C?( Pn - Pp )\ (13) 
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which trivially leads to the following relations between both sets of the coupling constants 



CL = C?-Cl (14a) 

CLm = 2Cf, (14b) 

or equivalent ly 

Co = 2^sum + ^tot) (15a) 

C{ = \C^ nm . (15b) 

Analogous relations hold for all other terms in the energy density. 

The Coulomb energy is a sum of the direct and exchange contributions, 



cCoul cCoul I cCoul / -i c\ 

° — c dir ' exch ' \ LU ) 



for 



fg- = e l f^ ftW^ , (17a) 

« = J J^dr/-^p^lA, (17b) 

where p p (r) and p p (rx, r 2 ) are the local and non-local proton densities, respectively, Eqs. 
(|5aD and (||), and e is the elementary charge. By using expressions (17) we approximate 
the charge density by the proton density without taking into account neither proton and 
neutron charge form factors [25| nor other corrections [26] which turn out to be small. 

In the algorithm, the direct Coulomb energy is calculated as the trace of the proton 
density matrix with the HO matrix elements of the Coulomb potential (see Sec. [5]). The 
Coulomb energy £7 Coul (r) is given by 

U c ™ l (r)=e 2 I dr'M^X, (18) 



where, in order to express ?7 Coul (r) in units of energy, the additional factor e is added as 
compared to the standard electrostatic expression. 



The exchange Coulomb energy is calculated in the Slater approximation [27, 28 



£^ = ~ 3 4& /3 ldr P i /3 (r). (19) 



4 \7T 

The term £ pa,n is equal to the average value of the seniority pairing interaction [JT| 
calculated in the BCS state: 

£ pair = - E (20) 

a=n,p \ i / 
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where G n and G„ are the neutron and proton pairing strengths, respectively, taken from 



Ref. |[29|| . Sums are performed separately over the neutron and proton single-particle 
states, while the standard BCS occupation factors, v^ a satisfy as always vf a + uf a =l. 
Note that £ pair does not correspond to the difference between the energies of unpaired 
and paired solutions, because the occupation factors also influence other terms in the total 
energy. In the present version of the code hfodd the pairing correlations are included 
only within the non-rotating case, see II. 

2.3 Constraints 

The total energy (||) can be minimized under specific constraints. One of the most essential 
ones is related to the so-called cranking approximation which is equivalent to solving the 
time-dependent HF equations in the laboratory frame or the HF equations in the rotating 
frame [p]]. To distinguish between the energy operators and/or eigenvalues written either 
in the laboratory frame or in the turning reference frame, the latter ones are referred to 
as Routhians. This notion applies also to situations when other constraints are taken into 
account. 

In order to find a constrained minimum of energy one has to find a minimum of the 
Routhian £' defined by 

£' = £ + £ mult + £ cran + £ numb , (21) 

i.e., equal to the sum of the energy and the terms responsible for the multipole, cranking, 
and particle-number constraints, respectively, as defined below. 



The multipole constraints are assumed in the standard quadratic form |30[ : 



£" mXt = Y.Cx»((Qx»)-Qx») , (22) 

A/i 

where {Q\^} are the average values of the mass-multipole-moment operators, Qx^ are the 
constraint values of the multipole moments, and Ca m are the stiffness constants. 

The cranking constraints aim at finding solutions corresponding to non-zero angular 
momenta, and are assumed in a combined form of the linear and quadratic constraints: 

£ CTan = -LUj(j y ) + Cj ((Jy) -Jy)\ (23) 

where J y is the operator of the component of the total angular momentum along the y 
axis, and J y is the corresponding target value. The code hfodd uses the y axis as the 
cranking axis and assumes that the x-z plane as the conserved symmetry plane (see below 
for more details). For the combined constraint, the physical angular frequency u is given 
by 

0{Jy) 

It turns out that a pure linear constraint (Cj=0) leads to a much more stable convergence 
properties, and then uj y =ujj. 
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The particle-number constraints ensure the correct average values of neutron and 
proton numbers when the pairing option is used, and are assumed in the standard linear 
forms, i.e., 

e™ a *> = -\ n (N n )-X p (N p ), (25) 

where A n and X p are the neutron and proton Fermi energies, and N n and N p are the 
corresponding neutron and proton particle-number operators. The particle number con- 
straints are taken into account only when the pairing correlations are included (in the 
non- rotating case, at present). Otherwise the numbers of particles are defined by select- 
ing a given number of occupied states explicitly (see II, Sec. 3.4). 



2.4 Hartree-Fock mean fields 



Upon a variation of the total Routhian (21) with respect to neutron and proton single- 
particle wave functions, one obtains the HF single-particle Routhian operators in the 
form: 

ft 2 

K = -— a + (r c vcn + r° dd + rf en + r° dd ) + u mult - u y j y , (26a) 
h' p = -—A + (r* ven + r° dd - rf 611 - r° dd ) + u Coul + u mult - u v j y . (26b) 



The momentum dependent potentials V are given by [[22|, 2^| 

1 

2/ 



r even = _v • [M t (r) V] + t/ t (r) + ^( W • S t (r)+ S t (r)- W ) , (27a) 



odd 



cr-CAr) V 



<x ■ S t (r) + i(V ■ I t (r) + I t (r) ■ v). (27b) 



2i 

In Eqs. (26), £/ mult represents the terms originating from the multipole constraints: 

[/-ult = 2 £ ^ ( { q x J _ (2g) 

With pairing correlations taken into account, the single-particle Routhians should in prin- 
ciple also contain the terms —X a N a originating from the particle number constraint ([25]). 
However, usually these terms are kept aside and included in the equations explicitly. Such 
a convention allows keeping for the single-particle Routhians the convenient standard en- 
ergy scale. 

The functions defining the HF potentials (27) depend on the local densities (5)-(7), 
and read [Q 

U t = 2C t p Pt + 2C t Ap A Pt + CjT t + C7 J V-J t + Ul, (29a) 
S t = 2C°s t + 2Cf s As t + CfT t + c7'Vxj t , (29b) 
M t = Clp t , (29c) 
C t = C?8 t , (29d) 

B t = 2Ctj t -C? J VPt, (29e) 
I t = 2CU t + C7 j V x s t . (29f) 
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The tensor gradient operators in Eqs. (27) and Q29ej ) are denned |22| as (Vcr) A1!/ =V At cr I/ 
and V^=£ A e^ A V A . 

In Eq. ( |29a| ) the term U' t represents the rearrangement terms resulting from the density 
dependence of the coupling constants. In standard parametrizations of the Skyrme forces, 
which are implemented in the code hfodd, only the Cf and Cf coupling constants depend 
on the isoscalar density p , and then 



u't = s* E iF-d + ^4 • (30) 




2.5 Hartree-Fock equations 

The eigenequations for the HF single-particle Routhians (26) are called HF equations, 
and read 

ti a ipi, a (ra) = e' ia ij) ija {ro). (31) 

The single-particle index i comprises all nucleonic quantum numbers. We use the conven- 
tion that this index has different values for both members within the time-reversed, or 
signature-reversed, or simplex- reversed couples (see Sec. i.e., every single-particle state 
is considered separately, and the sums over i are performed over all neutron or proton 
single-particle states. This is in contrast to some formulations where the sums are per- 
formed over only one-half of states which have one common value of, e.g., the signature, 
or one sign of the spin projection. [Note, for example, the factor \ in Eq. (|20|) .1 

The standard way of assessing the quality of convergence of the HF equations is to 
compare the value of energy £ calculated from Eqs. (^)-(|20|) with that calculated by using 
the sum of the single-particle energies 

£ S ' P - = £<Aa. (32) 

A connection between these two energies exists because £ s p ' is related to the average value 
of the single-particle Routhian, which for a density independent interaction and with no 
constraints is trivially equal to the average kinetic energy plus twice the average two- 
body interaction energy [jlj. In the realistic case of the density-dependent interactions, 
this simple relation has to be modified and one defines the equivalent energy £ by 

C J_$?s.p. I ^lckin I cpair crear _i_ ^cCoul crnult ccran ('qq'i 

° ~ 2 ' 2 ' ~ ' S^exch C corr — ^corr \ M ) 

Here, the rearrangement energy £ rear results from the density dependence of the Skyrme 
interaction and is equal to one half of the average value of rearrangement potential Eq. 

£ rcar = | / d 3 rU^ Po . (34) 

Similarly, the Coulomb exchange energy ( [L9f) can be considered as resulting from a zero- 
order interaction term depending on the density as pp 2 ^ 3 - Hence the resulting rearrange- 
ment term in Eq. (^) is equal to the Coulomb exchange energy with the factor |. Finally, 
the corrections resulting from the constraints read 

Crf = E ((Q*) - Qx») (Q*) (35) 
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for the multipole constraints (p8|), and 

^cott — -^yiJy) (36) 

for the cranking constraint with u y given by Eq. (^4j). Since the term corresponding to 
the particle number constraint is not included in the Routhian, the correction for this 
constraint does not appear in £. 

The difference between the energies £ and £ is exactly equal to zero when the densities 
and fields do not change from one iteration to the next one. Hence their difference 

5£ = £-£ (37) 

provides a useful measure of the quality of convergence, and is called the stability of the 
HF energy. In practice, £ approaches faster the final value of the total energy energy than 
does £. 



2.6 Hartree-Fock densities 

In terms of the Routhian eigenfunctions, the non-local density matrix is given by 

p a (ra,r'a') = E^H^JrV). (38) 

i 

The occupation factors vf a are determined through the standard BCS procedure (when 
the pairing correlations are taken into account for the non-rotating case of the code 
hfodd), or are equal to 1 or according to the chosen configuration (when the pairing 
correlations are not used). 

The energy densities (12) depend on the six real densities p t , r t , s t , T t , j t , and 

J t . Each of them has the isoscalar (t=0) and isovector (t=l) form, or equivalently, the 
neutron and proton form. Moreover, in Eqs. (12) and (29) there also appear six secondary 
densities which are some particular derivatives of the six basic ones, i.e., Ap t , V ■ Jt, As t , 
Vpt, V x j t , and V x s t . In principle, they can be calculated by an explicit numerical 
derivation of the six basic densities. However, the organization of the code hfodd allows 
calculating them directly at a negligible expense of the CPU time. This is so, because 
the principal numerical effort must anyhow be devoted to summing up the contributions 
of the HO basis states to the HF eigenstates at every point in space (see Sec. |42|) . Once 
this is accomplished, the remaining operations are very rapid. 

The scheme used in the code hfodd is therefore the following. First, at every point of 
space, r, and for every spin value a, the values of the HF wave functions, their derivatives 
in three directions, and their Laplasians are calculated. At this point it is convenient to 
combine the value of each wave function and its three derivatives into one four- vector. To 
this effect we introduce the symbolic notation for derivatives, V/i, which are distinguished 
by hatted indices taking values £i=0, 1, 2, 4. Then, /2=0 corresponds to the unity operator 
(no derivative involved), and the three usual indices correspond to the usual gradient in 
space, i.e., V^=(l, V). By using this notation we define the two following generic matrices 
of densities: 

D%, a = Y,< a [^^A^)]Wu^Ura')}, (39a) 
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L q / = £<Aa(™)A^* a (rV), (39b) 

i 

where indices q and q' denote the signs of spin variables a and a', respectively, i.e., q=+ 
for o— J r\ and q=— for a——\. 

The local densities required for the calculation of energies and fields have to be de- 
termined separately for neutrons and protons. Below we omit for simplicity the isospin 
indices and denote by 3? and 9* the real and imaginary parts. Then the required 34 real 
densities are given in terms of 64 complex or real densities Dpi and 4 complex densities 



L qq ' by the following expressions: 

• scalar densities 

p = D+ + + Dw, (4o; 

* = EK;+v). ( 4i : 

Ap = 23? [L ++ + L~) + 2r, (42; 

V-J = -29 (£> 2 + 3 - - £>+" ) - 25? (D+f - £>+") + 23 - , (43; 

• vector densities 

si = 23?D+~, (44a 

s 2 = (44b; 

s 3 = L^V-LW, (44c 

7\ = 23? ]>>+", (45a 

T 2 = -23 (45V 

T 3 = Efc + -^7). ( 45c 

Asi = 23? (L + - + L"+) + 27\, (46a 

As 2 = 23 (L+- - L _+ ) + 2T 2 , (46b; 

As 3 = 2 - L~) + 2T 3 , (46c; 

Vp M = 2St(D++ + D^), (47; 

= %(d; + + d; ), (48; 

(Vx«)x = 2» (£>++- Di") + 29f (£>+" + £>+-) , (49a 

(Vxs) 2 = 2» (£>+" + £>+-)- 2» (£>J^ - £>ro") , W 

(Vxs) 3 = -23(d+-+D V)-23?( j D 2 V + j D V) ) (49c 
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V x 


3 


h = 


23 




+" D 32~ 


V x 


3 


>2 = 


-23 




+■ 


V x 


3 


>3 = 


23 




+■ 



tensor density 



J, 



3 



J, 



^2 



J, 



^3 

3 Self-consistent symmetries 
3.1 Simplex 



&(d' 



D, 



J ' 



(50a) 
(50b) 
(50c) 

(51a) 
(51b) 
(51c) 



The reader is referred to Ref. ]3T[] for a detailed discussion of various possible choices of 
conserved symmetries and of the implications thereof. In the present version of the code 
hfodd we assume x-z plane as a priori the only symmetry plane of the HF states. Such 
an assumption has been dictated exclusively by the fact that the CPU time for a typical 
run is markedly lower when symmetries are assumed. However, in physical applications 
one usually requires that all symmetries are determined by the physical system itself as 
a result of the self-consistency condition and not by a pre-supposition. We find that the 
assumption of a single symmetry plane provides an acceptable compromise between the 
above conflicting criteria. 

As a consequence of the above assumption, the y-simplex operator, defined by 



S y = Pexp (— in X 



y i ' 



(52) 



where P is the parity operator, commutes with the single-particle Routhian operators 
(26), i.e., 

0, (53) 



k nJ Sy 



h' <? 



and that the many-body wave function is an eigenstate of S y : 



(54) 



There are some technical advantages of choosing of the ^/-simplex as the conserved sym- 
metry. One of them is related to the convenient phase relations which are discussed below; 
another one consists in the fact that the average values of multipole moments: 



are real, i.e., 



Qx, = (Qxp) = (*|QamI*> 

Q\fi = Qxij.- 



(55) 
(56) 
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This property results from combining the standard time-invariance condition [53], 

Ql» = (-I)^a,-,, (57) 
with the transformation law under the operation of the y-simplex: 

SlQxA = {-lYQx,-„- (58) 

Equations (|BT)|)-(|57]) allow to consider only (real) multipole moments with non-negative 
magnetic components, fi > 0. Apart from a restriction to real values, all multipole 
moments can be non-zero. In particular, deformations corresponding to four magnetic 
components, fi=0, 1, 2, and 3, of the octupole moment can be simultaneously present. 

The choice of the ^/-simplex as a conserved symmetry requires that the nuclear rotation 
(cranking) takes place around the y axis. In the following we refer to the y-simplex as 
the simplex, and keep in mind that this axis corresponds to the rotation axis. 

Because 

S 2 y = (-1) A , (59) 

the simplex has four eigenvalues: S—H and S=±i for even and odd numbers of particles, 
respectively. Due to the commutation relations 053]), the single-particle Routhians do not 
have non-zero matrix elements between states belonging to different eigenvalues of S y . 
The equations for the single-particle eigenvalues can, therefore, be solved separately for 
both values of the single-particle simplex s=±i, and the single-particle wave functions 
ipi,a(r(j) are labeled by these values of the simplex, i.e., 

S y ipi, a (r<7) = Siipi^ra). (60) 

In the above relation we do not introduce Sj as an extra label of the wave function since, 
according to our convention, all quantum numbers characterizing the single-particle states 
are contained within the index i. Throughout the code, input data, and output files we 
associate the sign "+" with the eigenvalue s=+i of the simplex operator S y , and the sign 
"— " with the eigenvalue s——i. 



3.2 Time-reversal 

The code hfodd solves the HF equations with or without the time-reversal symmetry im- 
posed, depending on the value of the input parameter IROTAT (see II, Sec. 3.3). However, 
the antiunitary time-reversal operator T 

f = -ia y K , f 2 = (-l) A , (61) 

where K is the complex conjugation in coordinate space, is used throughout to establish 
the phases of the HO-basis simplex eigenstates, Sec. |4~T| . 

Since the time-reversal operator commutes with the simplex operator, 

T Sy = SyT (62) 

the time-reversed single-particle states belong to opposite simplex eigenvalues. (Recall 
that the antilinear time-reversal operator changes the sign of the imaginary simplex eigen- 
value.) Therefore, only matrix elements for one of the eigenvalues, say s=+i, have to be 
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evaluated, while those for s=—i can be obtained by the time-reversal. One has only to 
keep track of the time-even and time-odd terms in the Routhians. When both kinds of 
terms are non-zero, as is for the case of rotation, the Routhians for both s eigenvalues have 
to be diagonalized. However, when the time-reversal symmetry is present (no rotation), 
only one diagonalization is sufficient. 



3.3 Parity and signature 

The code hfodd solves the HF equations with or without the parity symmetry imposed, 
depending on the value of the parameter ISIGNY chosen in the input data (see II, Sec. 
3.3). Since the simplex S y , Eq. (p2[) , is in the present implementation always a conserved 
symmetry, the case of conserved parity is in fact realized by requiring that the y-signature 
(signature for short), defined by 

R y = exp (—inJy) , (63) 

is a conserved symmetry. Then, the HF single-particle states are in addition labeled by 
the signature eigenvalues, r=±i, 

Ryipi, a (rcr) =riif) ija (ra), (64) 

while the simplex is a product of parity tt and signature r, s=nr, or equivalently the parity 
is a product of simplex and signature, n=sr*. As before, the signature quantum number 
Ti is not attached explicitly as a label of the wave function. The s=+i block splits into two 
parity-signature blocks (7r,r) = (+l,+i) and (— 1,— i), and similarly the s=—i block splits 
into (— and These four blocks can be diagonalized separately. Throughout 

the code, input data, and output files we associate the sign "+" with the eigenvalue r=+i 
of the signature operator R y , and the sign "— " with the eigenvalue r=—i. 

In the case of conserved parity and signature, the multipole moments Qx^ vanish for 
odd values of A, which is a consequence of the transformation law for multipole moments: 

P ] Qx»P = (-1) A Qa m - (65) 



3.4 T-simplexes 

In addition to the two important cases of either the simplex S y alone appearing as a 
conserved quantity, or of the parity and signature R y being conserved, one may have two 
other cases of conserved symmetries. These symmetries cannot give any new quantum 
numbers simultaneously with the simplex s, because none of the other simplexes, S x and 
S z , or signatures, R x or R z , commutes with S y . However, we may still have conserved 
antilinear operators, which do not provide quantum numbers, but restrict spatial proper- 
ties of solutions. The code hfodd can treat two such a symmetries [3T|, the x-simplex T 
and the z-simplex T , 



£J = TP exp (-mi,) , (Sj) =1, (66a) 
Sj = fPexp(-tiTj z ) , (Sj) 2 = l. (66b) 
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These symmetries are requested by activating in the input data (see II, Sec. 3.3) the 
options ISIMTX and ISIMTZ, respectively. The product of these two symmetries equals to 
the ^-signature, 

S T J T Z = Ry. (67) 

Therefore, if the Hamiltonian commutes with both of them, the signature (and hence 
parity) must also be conserved. 

The two T-simplexes transform the multipole moments in the following way: 

(5j)tQ v £j = (-1)"Q^, (68a) 
(SjyQ x ,Sj = (-l) x Qx,-„ (68b) 



which combined with the time-reversal and y-simplex transformation laws (pf\ ) and ([58]) 
gives, respectively, the following properties of multipole moments: 

Q\[i = {-ITQ^, (69a) 
Qx» = (-1) A+m Qa m . (69b) 

Therefore, commutation of S% with the Hamiltonian enforces the vanishing of all odd 
magnetic components of all multipole moments. Nevertheless, with broken parity, the 
even magnetic components of odd multipoles, like Q30 and Q32, can be nonzero. On 
the other hand, conservation of enforces the vanishing of odd magnetic components 
of even multipoles, and the vanishing of even magnetic components of odd multipoles. 
Therefore, with broken parity, the moments Q31 and Q33 can be nonzero. Of course, 
conservation of both T-simplexes leaves nonzero only even magnetic components of even 
multipoles. Table [I] summarizes the properties of multipole moments under different 
symmetry conditions. One should note that the broken parity may lead to the system 
which has the center of mass shifted from the origin, Qiot^O and/or Qny^O, while for 
the both T-simplexes broken it may in addition lead to a shape which is not in the 
principal axes of the quadrupole tensor, Q2iy^0. In order to avoid physically meaningless 
"deformations" which in fact correspond to a shift or to a rotation of the whole system, 



one has to use suitable constraints, see Sec. IO| and II - Sec. 3.7, to enforce conditions 



Qio=Qn=Q2i=0 in cases where the non-zero values are, in principle, allowed by broken 
symmetries. 

Condition (|67|) means that there exist a possibility [[H]] of breaking both antilinear 
symmetries, S'J and Sj, and still conserving the signature (and hence parity). Such a 
possibility is, however, excluded from the present version of the code. This is so because, 
it would have corresponded to a state with vanishing octupole moments, and nevertheless 
non-zero moments Q41 and Q43 present in the principal-axes reference frame (Q2i=0). 
The shapes of that type are very exotic indeed, cf. Ref. |33[| , and probably should be 



studied after an investigation of those corresponding to non-zero moments Q41 and Q43 
induced by the non-zero octupole deformation. The latter possibility is allowed in the 
code HFODD. 

4 Cartesian harmonic oscillator basis 

The code hfodd solves the HF equations by expanding the single-particle wave functions 
ipi{ra) onto the deformed HO wave functions ipn x n y n z ,sX ra ) m the Cartesian coordinates, 
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N x N y Nz 

Mr*) = E E E E C A WK (70) 

n x =0 n y =0 n z =0 Sz= _ 1_ 1 

Here A^ x , A^, and iV 2 are the maximum numbers of the HO quanta corresponding to the 
three Cartesian directions. However, as discussed in Sec. 2 of II, the sums over n x , n y , 
and n z are performed over the grid of points which form a pyramid rather than a cube. 
The HO wave functions have the standard form 

^n x n y n z ,s z iro) = lp nx (x)^ ny {y)lp nz (z)8 Sga , (71) 

where 

VwW = ^(e>-^, (72) 

and £,n=b dimensionless variables scaled by the oscillator constants 



= ^mwjh. (73) 
Polynomials H^(£) are proportional to the standard Hermite orthogonal polynomials 

^°)(0 = (v^2^!ptf„(0, (74) 



and normalized as ^ 

/ dZHV>(t)H§ ) {$e-* = 6 m .. (75) 

J — oo 

When convenient, we also use the standard bra-ket notation: 



\n x n y n z ,s z ) =ip nxnynziSz (ra). (76) 



4.1 Simplex basis 



The y-simplex symmetry operator ( 52" ) transforms the HO states ( [75] ) in the following 

way 

S y \n x n y n z ,s z ) = {-l) ny+ ^~ Sz \n x n y n z ,-s z ). (77) 

Since in the present implementation of the code hfodd, the simplex symmetry is always 
assumed, it is convenient to use the HO basis composed of states which belong to a given 
simplex, i.e., 

\n x n y n z , s=+i) = (i ny \n x n y n z , \) - r ny+1 \n x n y n z , -|)) , (78a) 
\n x n y n z , s=-i) = (-i Uv+1 \n x n y n z , \) + i~ ny \n x n y n z , -|)) , (78b) 

for which 

S y \n x n y n z ,s=±i) = (±i)\n x n y n z , s=±i). (79) 

Since the HO wave functions are real, the time-reversal operator ( |BT]) transforms them 
in the following way 

f\n x n y n z ,s z ) = (-1) \n x n y n z , -s z ) (80) 
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The relative phases of states ( 78a ) and ( 78b ) have been chosen in such a way that the 
time reversal simply flips the simplex: 



T\n x n y n z , s=±i) 



±\n x n y n z ,s= =F«). 



Having the relative phases established, we may still arbitrarily chose the absolute 
phases of, say, the s=+i simplex eigenstates. The choice in Eq. ( |78a| ) is made by consid- 
ering the antiunitary operator K 



K = Tier, 



1. 



(82) 



This operator does not act on the space coordinates and therefore conserves quantum 
numbers n x n y n z . Since it is an antilinear operator with the square equal to one, the 
phases in the spin space can always be chosen |35[] in such a way that all the basis states 
are its eigenstates with the eigenvalues being equal to 1. Since K commutes with T, such 
a choice of phase convention made in ( |78a| ) applies in fact to both simplexes, i.e., 



K\n x n y n z , s=±i) 



\n x n y n z ,s=±i). 



One should stress that K is not a conserved symmetry, and therefore the HF single-particle 
states do not have any particular symmetry with respect to this operator. 



4.2 Matrix elements in the simplex basis 

The properties of the operator K and the resulting phase properties allow to simplify the 
form the single-particle potentials (27) when evaluated between the simplex eigenstates. 
Indeed, for such a choice various terms in these potentials have matrix elements which 
are either real or imaginary, but not complex. As can be easily checked, in the time-even 
potentials the spin-independent terms and the terms proportional to a z are real, while 
the terms proportional to a x and a y are imaginary, whereas, in the time- odd potentials 
the spin-independent terms and the terms proportional to a z are imaginary, while the 
terms proportional to a x and a y are real. As already mentioned, the operator K is not a 
conserved symmetry of the Routhians, the HF wave functions do not have any specific K 
symmetry and are complex. 

Let O denote an arbitrary operator acting only in the space coordinates (independent 
of spin) and let \n x n y n z ) denote pure HO states (no spin component). Then we may 
summarize the properties of matrix elements in the simplex basis in the following way: 



(n x n y n z , s=+i\Oa \n' x n y n z , s=+i) = 

(n x n y n z , s=+i\Oa x \n' x n' y n' z , s=+i) = 

(n x n y n z , s=+i\Oa y \n' x n' y n' z , s=+i) = 

(n x n y n z , s=+i\Oa z \n' x n' y n' z , s=+i) = 

where the factors F{n) are defined as 

F(n) = { 



{n x n y n z \6\n' x n' y n' z ) F(n y - n' y ), (84a) 

{n x n y n z \0\n' x n' y n' z }F(n y + n' y + 1), (84b) 

(n x n y n z \0\n' x n' y n' z )F(n y +n' y + 2), (84c) 

(n x n y n z \d\n' x n' y n' z )iF(n y - n' y + 1), (84d) 

for n even, . , 
for n odd. 
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These properties allow omitting the calculation of matrix elements for which F=0. As 
mentioned in Sec. |3.2| , the matrix elements corresponding to s=—i need not be explicitly 
calculated. 



4.3 Matrix elements in the coordinate space 

For the Skyrme interaction, the Routhian operators (26) are differential operators com- 
posed of terms up to the second-order derivatives. Therefore, the spatial matrix elements 
of the type (n x n y n z \0\n' x n' y n' z ) can be calculated by using differential formulas for the HO 
wave functions (|72"D, i.e., 

= blHV>{tp)e-&, (86a) 
= 4H%(Qe-^, (86b) 

ft 

where #^(0 and ^ 2) (0 

are the (n+1)- and (n+2)-order polynomials of £ defined by 



= 2ni^ 1 (£)-£#f (£), (87a) 
#f(0 = (e 2 -2n-l)i?f (0- (87b) 

Matrix elements of all types of differential operators (up to the second order) can therefore 
be expressed through the expansion coefficients C* n ,(dd'), 

n+n'+d+d 1 

HW{i)H«\t) = £ C k nn ,(dd')Hl°\0, (88) 

k=0 

with 0<d+d'<2, and through the integrals involving the polynomials //^ (£) only. 
Coefficients C^ n ,(dd') can be in principle calculated recursively or explicitly [T5|], however, 
a method which is simpler and less prone to programming errors has been used in the 
code hfodd, namely, they are calculated numerically with the machine precision by 
employing the orthogonality relations (|75|) of Hermite polynomials and the Gauss- Hermite 
quadratures. From now on, we consider explicitly only those terms which do not contain 
derivatives, while all derivative terms can be treated analogously by using conditions ([88]) 
for d,<f >0. 

Concentrating on terms without derivatives, the matrix elements of an arbitrary func- 
tion 0(x,y,z) can be calculated as 

(n x n y n z \6\n' x n' y n' z ) = £ C£ n , (00) £ C n %, (00) £ C£ ni (00)O w „ (89) 

kx ky kz 

where 

O kxkykz = J d^ x d^ z O (& jj, %)Hg&)Hg(Z y )Hg&)e-e-G-Z. (90) 
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Integrals ( pO|) are calculated by the standard Gauss- Hermite quadratures []3"4" |. These 
integrals can be calculated exactly by noticing that for the Skyrme interaction the func- 
tions 0(x,y,z) are linear combinations [see Eqs. (29)] of densities (39), which in turn 
are quadratic in the HO wave functions and their derivatives. Therefore, alii integrated 
functions 0(x,y,z) have a form of polynomials W(x,y,z) multiplied by the typical HO 
Gaussian factors, i.e., 

O (x, y,z) = W (x, y, z)e~^ x ~^~^, (91) 



for ^=b^x^, cf. Eq. (|7^). By including the Gaussian factors explicitly in the Gauss- 
Hermite quadratures one can calculate the integrals of the remaining polynomials exactly, 
i.e., 



Ok 



J d&dtyd&W (ft, f , fc)^ ) (^)< ) (^)< ) (6)e-^- 2 ^- 2 ^ 



L 



x 'Hy ^ly 



Er lx \" r ly r lz w ( ^ m v m * i 



In the final formula, rji are the standard Gauss-Hermite nodes 
weights wi are included in the matrices G l k , i.e., 



(92) 

and the corresponding 
(93) 



The orders L x , L y , and L z of the Gauss-Hermite quadratures can be estimated from 
the maximum numbers of HO quanta included in the basis, Eq. ([70]) . These quadratures 
give exact results, also for terms depending on derivatives, for 



L, 



2N^ + 2, 



(94) 



or larger. 

It turns out that the CPU time required to perform the summations in Eq. ( |9"2"D is very 
small compared to other parts of the code. It is so provided the order of summations is 
as indicated in the formula, whence at most four-fold nested loops are required. A much 
larger effort is required to perform the summations in Eq. (|89|), where, even with the 
indicated order of sums, the seven-fold nested loops are necessary. Even if, in practice, 
the dimensions of these loops are rather small, this requires a substantial numerical effort. 
This is also the part of the code which is least amenable to vectorization. Such consider- 
ations indicate that the use of the Gauss-Hermite orders fl9~4"| ) ensuring exact integration 
is fully justified, even if these can be quite large for large numbers of HO quanta, and are 
certainly much larger than those usually used by other authors and codes. 

From the presented analysis it is clear that the summations of the HF densities (39) 
have to be performed only at the spatial points which are defined by the Gauss-Hermite 
integration. Moreover, the same derivatives of HO wave functions as in Eqs. (86) can be 



3 Except from the terms resulting from the density dependence 
terms (|l9|). The direct Coulomb terms are treated separately in Sec 



3C7| 



, and from the exchange Coulomb 
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used to calculate the derivative densities in Eq. (39). The spatial symmetries discussed in 
Sec. |] are employed to reduce the number of points in space where the densities have to 
be summed up, i.e., for three, two, or one symmetry planes only 1/8, 1/4, or, 1/2 points, 
respectively, are used. Table || summarizes the symmetries of different HF densities with 



respect to the y-z, x-z, or x-y plane |JIJ. These three symmetry planes are independent 
of one another and appear for the conserved symmetries S'J, S y , and Sj, respectively, 
although, as discussed above, not leading to the conserved quantum numbers (except if 
the parity is implied). 

5 Direct Coulomb potential 

Direct Coulomb potential is given by three-dimensional integral (|T8| ) and therefore differs 
from all the Skyrme mean-field potentials (29) which are simple linear combinations of 
densities. Various methods of calculating this potential are used in the existing HF codes. 
A direct integration of the Poisson equation can be performed in the one-dimensional cal- 
culations restricted to the spherical symmetry. In two-dimensional calculations performed 
for axially deformed nuclei such an integration is also possible [p6 |. In three dimensions 



one may either solve the Poisson equation by using the conjugate gradient method 137] or 



replace the Coulomb interaction by an integral over Gaussian interactions [|15| , |3q| . The 
latter method is very convenient for the Gogny interaction for which the matrix elements 
of Gaussians have to be anyhow calculated. 

In the code hfodd we have implemented the method based on using the Coulomb 
Green function [pj , which is particularly well suited for calculations employing the Carte- 



sian HO basis. This is so because the Coulomb potential ( fig) is then, in fact, not explicitly 
required; it is enough to calculate its HO-basis matrix elements. 



The direct Coulomb potential (18) can be expressed through the Dirichlet Green func- 



tion G£i(r, r') in the following way |BP] 

t/ Coul (r) = e 2 J y d*r>G D (r, r') Pp (r') - -L jf dh' ^^ U c °»\r'). (95) 

The first term is the volume integral over an arbitrary closed volume, while the second 
integral is performed over the surface of this volume. The Coulomb potential has to be 
known on the surface. The Dirichlet Green function fulfills the Poisson equation for a 
point charge and vanishes at the surface. In the surface term, the normal derivative is 
calculated with respect to the outward direction perpendicular to the surface. 

In the present application it is convenient to use in ([J3]) the volume in the form of the 
parallelepiped 

—D x < x < D x , 

-D y < y <D y , (96) 
-D z < z <D Z . 

Then, the required Green function can be expressed in a separable form as 

GD(r ' r)= w; 3 .£ iT^im ' (97) 
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where the functions / equal to sine or cosine ensuring the Dirichlet boundary conditions: 



, v f cos(J M x M ) for j M even, , 
A-VW 1 sin (J MXM ) f or ^odd, lJ? 



^ = ^^. 09) 



and 

In expression fl97|) , the sums over j x , j^, and j 2 are performed over all non- negative 
integer values. In practice, these sums have to be restricted to finite ranges. In the code 
hfodd we use the summation bounds 

0<] x ,] y ,] z <N Cou \ (100) 

and jV Coul can be chosen based on the following considerations: It is clear that when 
using the Green function for the calculation of the Coulomb potential one must effectively 
perform the Fourier transform of the proton density. Therefore, the corresponding wave 
vectors have to cover the region of momenta for which the proton momentum density 
distribution is large. This distribution is localized in the momentum space, i.e., does not 
extend to very large momenta; therefore, the sums can be terminated for indices giving 
suitable large fixed values of J M . The value of jV Coul depends, therefore, linearly on the 
steps in the momentum space given by 7r/(2D fl ), and hence inversely on the dimensions 
D^. Therefore, we are interested in working with appropriately small sizes of the 
parallelepiped. 

On the other hand, one has also to know the values of the Coulomb potential on the 
six faces of the parallelepiped. Suppose these faces are located far enough, where the 
proton density is already very small. Then the Coulomb potential on these faces can be 
very well approximated by the standard multipole expansion p9 



U C ° U \r) = £ (2A + 1)r 2A + i ^%,(6l, (10!) 

For the validity and precision of the multipole expansion truncated to a few lowest values 
of A, we are interested in working with appropriately large sizes of the parallelepiped. 

In practice it turns out that one can chose suitable dimensions of the parallelepiped so 
as to fulfill both requirements simultaneously. In the code hfodd the dimensions in three 
directions are chosen to be proportional to the corresponding oscillator lengths 1/6^, i.e., 



d_ 

V2K 



D » = ( 102 ) 



where d is a dimensionless parameter defining the overall size of the parallelepiped. Based 
on the performed numerical tests, the recommended values are d—20 and N ul =80, which 
for a heavy nucleus corresponds to the parallelepiped of about _D M ~30fm, and to the 
maximum momentum and the integration step in the Fourier transforms of about 4fm _1 
and 0.05 fm -1 , respectively. For the SD states in 152 Dy the maximum contributions of 
the monopole, quadrupole, and hexadecapole terms to the Coulomb potential are at the 
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transforms of the solid harmonics r^Y^ on the six faces of the parallelepiped. The values 
of proton multipole moments Q p x are anyhow routinely calculated for the determination 



surface of the D At ~30fm cube equal to about 3MeV, 50keV, and 1.5 keV, respectively. 
In the multipole expansion of Eq. (|101|) the multipoles above A=4 can therefore by safely 
neglected. 

The calculation of the HO matrix elements of the Coulomb potential can now be very 
easily completed. First, we insert the multipole expansion ( |101| ) in the surface term in Eq. 
J5|). All surface integrals can now be performed, because they amount to suitable Fourier 

of nuclear deformations. 

Second, we use the fact that the proton density has the form of a polynomial multiplied 
by the HO Gaussian factor [cf. Eq. ([91]) 1, 

(f (x, y, z) = W p p (x, y, ^e - ^-^-^, (103) 

and this polynomial can be expressed exactly as a linear combination of Hermite polyno- 
mials, 

W*(x,y,z)= £ ^U*.^^^)^^!/)^?^^)- (104) 

kx ky kz 

Here, for later convenience we have included the factors \/2 in the arguments of the 
normalized Hermite polynomials. 

Coefficients p' kxkykz can be calculated (again exactly) by using the Gauss-Hermite 
quadratures of the order given by Eq. (|94"D in the formula analogous to (Q). In fact, if we 
want to use the same Gauss-Hermite nodes as before, we have to employ the formula 

p'l k k = V G'i V G ,l l jr G'l* W* (jt, , (105) 

lx = 1 ^ y — 1 «Z = 1 

for 

G' l k = Wl Hf\ Vl ). (106) 

Third, we insert the proton density in the form of Eqs. ( [10 3D and ( [1041 ) in the volume 
integral in Eq. ( |95|) . Then we see that the Fourier transform of proton density amounts 
to calculating the Fourier transforms of the HO oscillator wave functions. In practice, 
the dimensions of the parallelepiped are large enough to replace these transforms by the 
analytic results pertaining to integrations extended to infinity, i.e., 



B{ = f" d V f(u 3 -r,)Hi 0) ( V )e 

J — d 



27r#£ 0) (c^)e-H(-l)[!] for k + j even, 

for k + j odd, 



(107) 



where 



= ^T^> ( 108 ) 



and 



(j + 1)tt 
2d 

denotes the integer part. Due to the choice of parallelepiped dimensions adapted 



to the basis, Eq. (|102|) , the factors B k are the same for each Cartesian direction. 
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After these manipulations, the Coulomb potential, the volume and surface terms alike, 
are expressed as linear combinations of terms proportional to f{J x x)f(J y y)f(J z z) . Instead 
of performing this sums, and calculating the Coulomb potential in spatial coordinates, 
one may rather directly calculate its HO-basis matrix elements. This gives 

q p 2 t/. . . _i_ y rf> q a m 

*■*>*• ndf k *f "yj^ «*bl(j x + l)z + bl(j y + iy + b*(j z + ir 1 1 
where the volume term reads 

Yjxjyj* = £ £ ^ky £ Bk z P k x kyk z l (HO) 

kx ky kz 

and the sums in Eq. (|109|) are restricted to the finite set of indices ( |100|) . 

In Eq. ( |109|) the surface term is a linear combination of proton multipole moments 
multiplied by fixed matrices ^ . These matrices are calculated only once by integrating 
solid harmonics on the faces of the parallelepiped. Finally, the Coulomb matrix elements 
are obtained from 

(n x n y n z \U c ™y x n' y n' z ) = £ £ <?* £ C^V^, (111) 



k z 



where the new set of coefficients C' k nn , is required to accommodate the additional factor 
of y/2, i.e., 

n+n' 

H$»(t)H$>{£) = £ C' k nn ,Hf\y/20. (112) 

fc=0 

Similarly as for other terms of Routhians, evaluation of sums in Eqs. (|109|) and (|110|) 
requires much less numerical effort than of those in Eq. ( |111| ). 

The Green function method allows calculating the matrix elements of the Coulomb 
field directly from the matrix elements of the proton density. The summations involved 
require the numerical effort which is typical for other terms of the Skyrme interaction, in 
spite of the fact that the Coulomb interaction is not of zero-order. Performing the sums 
in the order indicated in Eqs. ( |109| ) and (|110 ) never leads to more than four-fold nested 
loops over low-dimension indices. 
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Table 1: Pattern of allowed non-zero values (denoted by X) of the expectation values 
of the multipole moment operators for four different cases of symmetries treated by the 
hfodd code. 



Symmetry planes: 


one, x-z 


two, x-z and y-z 


two, x-z and x-y 


three 


Conserved: 


Sy 


S y and 5j 


S y and S^ 


Sy, «Sj, and Sj 


Qio 


X 


X 








Qn 


X 





X 





Q20 


X 


X 


X 


X 


Q21 


X 











Q22 


X 


X 


X 


X 


Q30 


X 


X 








Q31 


X 





X 





Q32 


X 


X 








Q33 


X 





X 





Qao 


X 


X 


X 


X 


Qai 


X 











Qa2 


X 


X 


X 


X 


Q43 


X 











Q44 


X 


X 


X 


X 
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Table 2: Symmetries of HF densities with respect to the three symmetry planes appearing 
when one of the three different symmetries are conserved. 



Conserved: 




Sy 


si 


Conserved: 




Sy 


S T Z 


Plane: 


y-z 


x-z 


x-y 


Plane: 




X-Z 


x-y 


P 


+ 


+ 


+ 


Ap 


+ 


+ 


+ 


T 


+ 


+ 


+ 


V J 


+ 


+ 


+ 


Si 






+ 


Vip 




+ 


+ 


S2 


+ 


+ 


+ 


V 2 p 


+ 




+ 


S2 


+ 






V 3 p 


+ 


+ 




?i 






+ 


Asi 






+ 


T 2 


+ 


+ 


+ 


As 2 


+ 


+ 


+ 


T 3 


+ 






As 3 


+ 






h 


+ 


+ 




(V x s)i 


+ 


+ 




h 








(V x a) 2 








h 




+ 


+ 


(V x a) 3 




+ 


+ 


Jn 








(V X J)! 






+ 


J21 


+ 


+ 




(V x j) 2 


+ 


+ 


+ 


J31 


+ 




+ 


(V x j) 3 


+ 






J\2 


+ 


+ 












J22 
















J32 




+ 


+ 










Jl3 


+ 




+ 










J23 




+ 


+ 










J 33 

















26 



